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Summary 

A line Gauss-Seidel (LGS) relaxation algorithm 
in conjunction with a one-parameter family of up- 
wind discretizations of the Euler equations in two 
dimensions is described. Convergence of the ba- 
sic algorithm to the steady state is quadratic for 
fully supersonic flows and is linear for other flows. 
This is in contrast to the block alternating direc- 
tion implicit methods (either central or upwind dif- 
ferenced) and the upwind biased relaxation schemes, 
all of which converge linearly, independent of the flow 
regime. Moreover, the algorithm presented herein 
is easily coupled with methods to detect regions of 
subsonic flow embedded in supersonic flow. This al- 
lows marching by lines in the supersonic regions, con- 
verging each line quadratically, and iterating in the 
subsonic regions, and yields a very efficient iteration 
strategy. Numerical results are presented for two- 
dimensional supersonic and transonic flows contain- 
ing oblique and normal shock waves which confirm 
the efficiency of the iteration strategy. 

Introduction 

A large class of problems of fundamental impor- 
tance to the field of computational aerodynamics is 
that of simulating the flow of an inviscid, compress- 
ible gas in the transonic and supersonic regimes. The 
most popular approaches used in obtaining steady- 
state solutions to these problems have been the im- 
plicit spatially split methods of the approximate fac- 
torization (AF) type (refs. 1 and 2) and the highly 
vectorizable explicit schemes such as the multistage 
Runge-Kutta method used by Jameson and Bakes 
(ref. 3) and the predictor-corrector scheme of Mac- 
Cormack (ref. 4). A motivation for algorithm re- 
search stems from the fact that the highly vector- 
izable explicit schemes suffer from stringent stabil- 
ity restrictions, and, though the implicit AF meth- 
ods are unconditionally stable in two dimensions, 
they require an optimal set of iteration parameters 
for rapid convergence which are difficult and time- 
consuming to obtain. (See ref. 5.) Moreover, the 
implicit central-differenced AF schemes are unstable 
for three-dimensional Euler equations. (See ref. 6.) 

The recent emergence of upwind-differencing 
technology for the Euler equations has opened the 
door to a new class of solution strategies aimed 
at improving computational efficiency. Since up- 
wind differencing is more costly per iteration to im- 
plement than central differencing, overall efficiency 
gains must be derived from improved convergence 
rates. This is made possible by the improved condi- 
tioning of the system of difference equations afforded 
by the more physical upwind discretization. 


Some preliminary work by Chakravarthy (ref. 7) 
and Van Leer and Mulder (ref. 8) using relaxation 
methods and upwind schemes has been promising. 
In the present study, the fact that upwind differ- 
encing closely models the propagation of information 
along characteristics is exploited by choosing a com- 
bination of line relaxation and upwind discretization 
such that, for fully supersonic flow in the stream- 
wise (marching) direction, the algorithm becomes a 
direct solver of the linearized problem. That is, the 
algorithms of the family presented here all revert to 
an efficient implementation of Newton’s method in 
the supersonic regime and result in quadratic con- 
vergence to the steady state. 

Further efficiency gains are realized in supersonic 
regions by iterating on each line to remove the lin- 
earization error before proceeding to the next line, 
thus solving the problem in a marching fashion. 
This can be accomplished since, in the absence of 
waves propagating adverse to the marching direction, 
the combinations of line relaxation and upwind dis- 
cretization presented in this study result in the un- 
coupling of the discrete algebraic problem by lines. 
Thus, the computational efficiency of space marching 
methods can be effectively recovered. For problems 
with mixed supersonic-subsonic regions, the solution 
strategy can be extended to detect the flow-regime 
interfaces, march in the supersonic regions, and iter- 
ate with line Gauss-Seidel (LGS) in the subsonic re- 
gions. A comparison of CPU times with and without 
the new iteration strategy is provided in the section 
entitled “Results.” 

The authors would like to acknowledge several 
fruitful discussions on upwind differencing with J. L. 
Thomas of Langley Research Center and Bram van 
Leer of the Delft University of Technology, in the 
Netherlands. 

Symbols 

A, B Jacobian matrices 

c speed of sound 

e total energy per unit volume 

/, g mass, momentum, and energy fluxes 

/ identity matrix 

y, k x,y grid-line indices 

M local Mach number 

p pressure 

q conserved variables 

R residual, Riemann invariant 

5 entropy 


t 

time 

U, V 

x, y Cartesian velocity components 

U 1, C/2 

upper diagonals 

x,y 

Cartesian coordinates 

1 

ratio of specific heats 

A, V 

finite-difference operators 

Ax, Ay 

mesh increments in x- and y- 
directions 

8q 

change in q over a time increment 

6t 

time increment 

K 

spatial differencing parameter 

V 

Mach number switch 


• Dependent-variable evaluation 
■ Flux evaluation 



p density 

4> spatial differencing switch 

Subscripts: 

A , B, L, R state quantities from above, below, 
left, and right, respectively 


In equation (2), 
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j,k spatial indices (1 to J and 1 to K) 

x, y x, y coordinate directions 

oo free stream 

Superscripts: 

± positive and negative flux 

contributions 

n iteration index 


Spatial Discretization 

The integral form of the conservation equations 
governing the flow of an inviscid fluid can be written 
as 

~ [ qdV + j F -nds — 0 (1) 

ot Jv Js 

where V is the volume bounded by the surface S, n is 
the outward-pointing unit normal to the surface, q is 
the vector of conserved-state variables, and F is the 
flux vector. For simplicity in describing the upwind 
method, consider the semidiscrete formulation of 
equation (1) for a uniform Cartesian grid. 

Sketch A shows the notation for the semidiscrete 
formulation. Dividing by the cell volume (Ax, Ay) 
and rearranging yields 



_ (/j+l/2,k -/j-l/2,fc) 

+ Ky ( 9 >’ k + 1 / 2 ~ 9 i< k ~ 1 / 2 ) 


~ R j,k 


( 2 ) 


where p is the density, u and v are the velocity 
components in the x- and y-directions, respectively, 
p is the pressure, and e is the total energy per unit 
volume. The equation set is closed by the ideal gas 
law. 

In this study, upwind differencing is implemented 
by splitting the convective fluxes, / and y, into pos- 
itive and negative contributions. The Jacobian ma- 
trices associated with the linearization of the positive 
split-flux contributions have eigenvalues greater than 
or equal to 0; likewise, the Jacobian matrices of the 
negative split-flux contributions have eigenvalues less 
than or equal to 0. The particular flux vector split- 
ting (FVS) technique used here is that of Van Leer. 
(See ref. 9.) Although other FVS techniques could 
be implemented in the general framework which fol- 
lows, Van Leer’s was chosen because the individual 
split-flux contributions vary smoothly across sonic 
and stagnation points. Also, an important reason 
for choosing Van Leer’s scheme is its simplicity. The 
splitting of / is given by 

/(<?) = / + (<?) + /“(<?) (4) 

where, for \M X \ < 1, 


/*(<?) = 


f± 

J mass 


/massK'Y — l) u i 2c]/')' 
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L/massdh - 1)« ± 2 c ] 2 /[ 2( 7 2 - l)] + V 2 /2}_ 
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and 


/mass ” P c 


2 (Mx ± 1) 


speed. In the next section, a marching method for 
supersonic regions is developed based on a particular 
choice of k x . 


For M x > 1, /+ = / and /~ = 0; for M x < -1, /+ = 
0 and /“ = /. 

In the preceding series of equations, M x is the 
local Mach number based on the velocity component 
in the ^-direction u and on the local speed of sound 
c. The splitting of g in terms of M y = v/c follows 
similarly. 

Utilizing the flux vector splitting technique, a 
pair of state vectors are assigned to a common cell 
interface, and a single numerical flux is derived from 
this pair. For example, fj+\/2,k can b e expressed as 

/(?)i+l/2,fc = / (QLM)j+l/2,k 

= f + ( QL)j+i/2,k + f~ (<}R)j+l/2,k ( 5 ) 

where the subscripts L and R imply that the de- 
pendent variables at the cell boundary are evaluated 
from the left and right, respectively. The estimates 
for the state vector pair at the (j + 1/2, k) cell bound- 
ary are given by 

(9l)j + 1/2,* = \ K 1 - **) V + (1 + «*) A]gy |fc (6) 

(QR)j+l/2,k = 9j'+l,fc ~ ^ K 1 - K x) A + (1 + *x) ^Ify + l,* (7) 

where 

= Qj +l,k ~~ Qj,k (®) 

^Qj,k = Qj,k ~ Qj—l,k (9) 

Similar formulas hold for the interpolation in the 
y-direction with parameter K y . For 0 = 0, the 
method is first-order accurate in space, and, for 
0 — 1, the parameters k x and K y control the spatial 
accuracy. This accounts for the truncation error of 
the method. Some common choices are k = 1/3, 
which corresponds to the only third-order accurate 
member of the family (strictly for a one-dimensional 
problem) and k = — 1, the fully upwind second-order 
scheme. The remaining second-order members are 
all upwind biased or centered approximations. 

To reduce oscillations near discontinuities, a 
limiter may be incorporated into the scheme by 
modifying the interpolations given in equations (6) 
and (7). (See ref. 10 for two examples.) In this study, 
a monotonicity-preserving limiter was not used, but 
the effect of limiting was emulated by reducing the 
second- and third-order interpolations of q on a cell 
face to first order if the higher order interpolation re- 
sulted in a negative value for the square of the sound 


Relaxation Algorithm 

Application of the Euler implicit time-integration 
scheme in delta form and time linearization to the 
system of semidiscrete conservation laws given by 
equation (2) yields 


where 


L — 

6t dq 


$qj,k — Rj,k 


(10) 


6q = q n+1 — q n ; 6t n = £ n+1 — t n ;n = Iteration index 


In matrix notation, equation (10) can be expressed 
as 

N6q = —R (11) 

where TV is a large, banded, block coefficient matrix 
with a block size of four. For 6t tending to infin- 
ity, the solution of equation (11) by direct inversion 
is Newton’s method. In general, ultimately because 
of the bandwidth of TV, the direct solution of equa- 
tion (11) is not efficient. For a fixed grid, the band- 
width of TV depends on the particular values of 0 and 
k used in evaluating the state vectors at a cell inter- 
face. Even more important, however, is the fact that 
for the upwind methods, the bandwidth is also de- 
pendent upon the flow regime of the particular prob- 
lem under consideration. This is in contrast to stan- 
dard central-differenced schemes, which maintain the 
same structure of the coefficient matrix independent 
of the Mach number. 

To be specific, equation (10) can be expanded to 
read 


(8) 3 , fc + i{[ A+ [qL)SqL 

+ A ~ (9 « )l59 «]j + l/2,fc 


where 


' [ A ^ (<1l) S < 1L + A teR) S Q&\j-i/2 t k J 

+ ^ { [ S+ Sq B + B ~ (9A) M iifc+1/2 

- [ B+ {ib) 6 9B + B ~ (,1a) 6< Ia] j } = ~ R j,k 

(12a) 
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B ± = 
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dq 

3g± 

dq 


(12b) 
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2,3 
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2, K 


J, K 


The subscripts B and A imply the evaluation 
of q at the cell boundary from below and above, 
respectively. The structure of N is depicted in 
sketch B, using equations (6) and (7) and similar 
formulas for the state- variable interpolation in the 
{/-direction. Any individual element on the nonzero 
diagonals represents a 4 X 4 matrix. Vertical-line 
Gauss-Seidel (VLGS) relaxation sweeping through 
lines in the positive x-direction would solve the linear 
system by replacing the coefficients on the upper 
diagonals, labeled U 1 and 172, with zero. With the 
ordering shown in sketch B, the elements on the 
upper diagonals are given by 


U 1 = 


1 

Ax 
+ A 


' A + ^ + A- [I 

3+1/2, k 4 + j'+l/2,k 1 1 

0(1 ~ Kx) 

j-l/2.k 4 


U2 — — — — 

Ax 


3 + 1/2,*' 


0(1 - «x)l 


_ <I>K X \ 

2 ; 

(13) 


(14) 


Supersonic Flow 

For a fully supersonic flow in the streamwise (x) 
direction A~ is zero everywhere; thus, 172 is zero. 
However, Ul is in general nonzero because of the 
presence of A + . By choosing 0 = 0 (first-order 
scheme) or 0 = 1 and k x = — 1 (fully upwind second- 
order scheme), the elements on Ul become zero, and 
vertical-line Gauss-Seidel becomes a direct solver of 
the linear problem, resulting in rapid convergence 
to the steady state. Even for subsonic flows, when 
the upper diagonals are nonzero, the line relaxation 
algorithm with alternating sweep directions results 
in an efficient solution method. (See refs. 7, 8, and 
ii.) 

For supersonic problems, the classical implemen- 
tation of VLGS or global iteration (solving the sys- 


tem of equations once on each line and then pro- 
ceeding to the next line) is not the most efficient 
approach. A superior implementation of VLGS is 
to start at the first interior line of the grid, next to 
the inflow boundary, solve the linearized problem, 
update the dependent variables and residuals, and 
continue the process until the machine zero steady- 
state solution on the first line is obtained. The next 
step is to proceed to the next line and continue until 
the last line is completed, thus recovering the steady- 
state solution. This strategy is referred to as local, 
rather than global, iteration. In this approach, as 
8t — ► oo, Newton’s method is being implemented by 
lines, so that the solution along each line converges 
quadratically. Note that the global strategy attempts 
to obtain the solution at a downstream field point be- 
fore a fully converged solution has been obtained up- 
stream. Based on the mathematical theory of char- 
acteristics, this is not possible, because the steady 
solution along the downstream lines depends on the 
steady state upstream. This dependence manifests 
itself in the fact that local iteration is more efficient 
than global iteration for this class of problems, al- 
though both converge quadratically. 

Transonic Flow 

For problems containing separate regions of su- 
personic and subsonic flow, a hybrid strategy com- 
bining local and global iterations can be developed. 
For this type problem, where large subsonic regions 
can exist, global iteration may be used with back- 
and-forth sweeping. A more efficient hybrid strategy 
can be developed if the steady-state flow-regime in- 
terfaces can be located. As an example, consider the 
dual-throat inlet configuration in sketch C, a pos- 
sible Mach number distribution for the inlet, and 
the coefficient matrix that would be obtained from 
a first-order upwind discretization of the problem, 
all of which are related to the iteration strategy that 
follows. 

Assume that the steady-state location of the 
two interfaces separating the supersonic-subsonic 
states are known. Since the upper diagonal in 
the coefficient matrix has nonzero elements only 
on the rows corresponding to the subsonic region 
and along the preceding supersonic-subsonic inter- 
face (but not vice versa), the problem can be solved 
by local iteration in the supersonic region until a 
supersonic-subsonic interface is reached. Then, since 
the beginning and ending rows belonging to the 
submatrix are known, the large submatrix depicted 
in figure 3(c) can be solved iteratively until conver- 
gence. The last row included is the subsonic line 
just preceding the supersonic-subsonic cell interface 
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at the second throat. It is sufficient to close the sub- 
matrix with this row, because there is no influence 
on the solution from the downstream supersonic line. 
Solving the submatrix in this fashion is referred to as 
block iteration, because it involves the iterative solu- 
tion of a group or block of lines in a global fashion. 
Having obtained the steady-state solution in this re- 
gion, the remaining supersonic field can be solved by 
local iteration. 

The automatic implementation of the hybrid 
strategy is very straightforward. All that is needed is 
a switch that determines whether a line can be solved 
locally or must be included in a block-iteration group. 
In sketch C, the parameter v serves this purpose, and, 
for the first-order scheme, is defined by 

f 1, if M(q R ) j+l /2,k > 1 ^ k (1 5) 

J \ 0, otherwise ' 



If Vj — 1, the line can be solved locally, otherwise 
it must be included in a larger submatrix. The value 
of Vj at the last line depends on the type of boundary 
condition specified at the outflow. Obviously, for a 
supersonic outflow boundary condition, Vj = 1; for 
subsonic outflow, Vj = 0. 

The situation is only slightly more complicated 
with a second-order, fully upwind approximation 
in the streamwise direction. The only additional 
requirement for local iteration on line j is that 
M( qR,)j-i/2,k — 1? this ens ures that the elements 
on the diagonal Ul in sketch B are zero. Therefore, 


the definition of Vj for <j> = 1 is 



if M{q R ) J+l/2 k > 1 

and — * f° r (16) 

otherwise 


In practice, the steady-state location of the 
supersonic-subsonic interfaces (and vice versa) are 
not known a priori, but must be determined as part 
of the solution process. Two approaches were tried 
in this study. The first approach was to apply the 
global iteration technique first until Vj became fixed 
in time and then to apply the line-block iteration 
technique for the final pass through the mesh. A 
slightly more sophisticated approach is to initially 
solve the problem on each line to a loose tolerance 
(e.g., a two- to three-order reduction in the lineariza- 
tion error and/or a fixed number of iterations). After 
one left-to-right sweep through the mesh, the values 
of Vj are computed and the line-block-iteration strat- 
egy is applied on the next pass. Subsequent hybrid 
iterations are made with the line-block strategy un- 
til Vj becomes independent of time. A machine zero 
tolerance is then set, and the steady-state solution 
is recovered on the next hybrid sweep through the 
grid. It is not known under what conditions (toler- 
ance level, problem severity, etc.) this latter strategy 
will succeed or fail. 


Initial Conditions for Local Iteration 

To start a problem, a simple set of initial condi- 
tions are usually specified (e.g., uniform flow at free- 
stream conditions everywhere). However, in super- 
sonic regions where local iteration is used, another 
choice that may reduce the computational effort is 
to use the solution from the previous line as an ini- 
tial condition. Although this choice may be an im- 
provement over that of the free-stream conditions, 
it is generally not a good enough guess to put one 
inside the domain of attraction of Newton’s method 
without requiring additional iteration. This situation 
can be at least partially alleviated with only a minor 
coding effort by using an initial condition generated 
by solving the steady-state Euler equations approxi- 
mately. As shown in the section entitled “Results,” 
significant savings are obtained by this approach. 

Implementation of this procedure is as follows. 
Consider the following differential form of the steady- 
state Euler equations in conservation-law form: 


d -l + d JL =0 

dx dy 


( 17 ) 
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Let 


Ajh = hj + 1 - hj ( h = /, g , or <?) (18) 

where the k subscript has been suppressed for 
simplicity. Using first-order differencing in the 
streamwise direction and using equation (17) yields 

v +a ^=- ai (I ) 3 (19) 

Linearization of A jf and A 3 g with respect to q yields 


df\ . d {dg 
dqji + ^dy (dg 


A ~q - - Ax 


dg\ 

dy) 


(20) 


In equation (20), the Jacobian matrices, and 


are already formed (eq. 12(a)); consequently, the 
implementation and solution of equation (20) is very 
straightforward. The spatial discretization in the y- 
direction is upwinded, as previously described. Then, 


9j + 1 - 9 ? + ( 21 ) 

is set as the initial condition on the j + 1 line. 


Operation Count, Time Stepping, and 
Boundary Conditions 


increase is certainly worth the effort if quadratic con- 
vergence, as opposed to linear convergence, can be 
obtained. 

The general time-stepping scheme is a variation 
of that described in reference 8. At any level of 
iteration, the time step is obtained from 

_ bt 0 

\\R n h 

where ||i? n ||2 is the L 2 -norm of the residual (the 
right-hand side of eq. (2)) normalized by the initial 
residual and St Q is the initial time step. A com- 
mon choice for 8t 0 would correspond to a maximum 
Courant number of the order of 10. This choice would 
allow the initial solution to adjust to the enforce- 
ment of the steady-state boundary conditions with- 
out incurring severe oscillations. As the residual ap- 
proaches zero, very large time steps are obtained to 
mimick a steady-state Newton method. 

Finally, it has been tacitly assumed that all 
boundary conditions are treated implicitly, which 
is required if quadratic convergence is to be ob- 
tained. Nonlinear, delta-form boundary conditions 
can be readily linearized by the techniques already 
presented and built into the discrete algebraic prob- 
lem to yield a completely consistent formulation of 
Newton’s method. 

Results 



Before proceeding to the test problems, a few ad- 
ditional remarks on overall implementation are in 
order. Frequently, questions arise concerning the 
computational effort associated with the solution of a 
system of equations characterized by a block pentadi- 
agonal coefficient matrix as opposed to a block tridi- 
agonal solution process. This is pertinent, since all 
the discretizations with 0 = 1 require a pentadiago- 
nal solver. The operation count, without pivoting, is 
approximately n[4m 2 + 2 m - 1] for a block tridi- 
agonal solver and n[9m 2 + 3m — 1] for a block 
pentadiagonal solver, where n is the number of un- 
knowns (mesh points) and m is the block size. (See 
ref. 12.) Thus, approximately 2.18 times more effort 
is required for that part of the iteration associated 
with the pentadiagonal solver than with a tridiagonal 
solver. The actual increase in effort has been mea- 
sured on a Control Data CY170-855 computer, and a 
2.154 ratio was obtained. Solving the linear system, 
however, represents only 30 percent of the computa- 
tional time per iteration. The remaining 70 percent 
is the effort required to set up the linear problem and 
the boundary conditions. The total increase in CPU 
time per iteration is therefore only 35 percent. This 


To assess the relative performance of the various 
iteration strategies, the same test problem is solved 
by different means and CPU times are compared. 
The first problem considered is the inviscid shock 
reflection problem (ref. 13), in which a shock wave 
of prescribed strength reflects off a flat plate. The 
problem is simulated by prescribing free-stream con- 
ditions at the inflow boundary (M 0 0 = 2.9) and 
extrapolating conditions at the supersonic outflow 
boundary (first order). Along the upper boundary, 
all conditions are fixed through the incident shock an- 
gle of 29°, which corresponds to an overspecification 
of the boundary conditions by setting the exact solu- 
tion. The boundary conditions on the plate were flow 
tangency and first-order extrapolation of />, u , and e. 
The grid contained 61 equally spaced points in the re- 
direction, 0 < x < 4.1, and 21 equally spaced points 
in the y-direction, 0 < y < 1. Pressure contours 
for the first-order x and y scheme (0 = 0), the fully 
upwind second-order x and y scheme (0=1, k x = 
— 1, K y = — 1), and the hybrid second-order x , third- 
order y scheme (0 = 1 , k x = — 1 , K y ~ 1/3) are 
compared in figure 1. Clearly, a dramatic resolution 
improvement is obtained with the 0=1 schemes. 
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The convergence rates of the three methods us- 
ing the global iteration procedure are depicted in 
figure 2. This problem was then solved by lines using 
the initial conditions given in table 1. 

Table 1. Comparison of CPU Time for 
Shock Reflection Problem 1 

[Tolerance = 10 1 3 ; L 2 -norm of residual] 




CPU time, sec, for — 



Global 

Uniform 

Previous 

Eq. (20) 

Eq. (20) 


Ky 

iteration 

flow 

line 

«» = 0) 

(0=1) 

0 


20.6 

7.9 

7.8 

3.8 

3.8 

1 

-1 

30.8 

(2) 

10.6 

6.5 

(2) 

1 

1/3 

40.8 

(2) 

11.8 

11.9 

12.4 


1 Control Data VPS 32 Computer; FORTRAN 2.1.5; OPT = 
BE compilation; scalar code. 

2 Did not converge because of state vector interpolation near 
shock reflection. 

The results shown in table 1 verify that local it- 
eration is significantly better than global iteration 
for supersonic flow. Using free-stream initial condi- 
tions on a line does lead to difficulty with the higher 
order methods if discontinuities are present because 
of the nature of the interpolation. Using the solu- 
tion from the previous line as the initial guess alle- 
viates this problem, but this method can be further 
improved by using the first-order-accurate solution 
of equation (20) as initial conditions. Finally, for 
the higher order methods, one might expect that us- 
ing a consistent y-discretization in equation (20) to 
generate the initial condition would be superior to 
the first-order discretization, but this has not been 
the case. The third-order discretization resulted in a 
marginal CPU increase, and the fully upwind second- 
order y-discretization resulted in failure. This failure 
to achieve a further CPU reduction was probably due 
to the fact that equation (20) is only first-order ac- 
curate in x. 

The next problem considered is the supersonic 
and transonic flow in a dual-throat inlet. The nozzle 
contour is very mild and was generated by specify- 
ing a parabolic Mach number distribution between 
the first throat and the exit and by then solving for 
the area distribution from the one-dimensional area 
Mach number relation. Ahead of the first throat the 
area was constant. A 69 x 31 mesh was used with 
equal spacing in the x-direction and slight cluster- 
ing in the y-direction. The upper-wall contour is 
sketched in figure 3 along with the converged center- 
line Mach number distribution from the first-order 
scheme for both supersonic (Moo = 1.10) and tran- 
sonic (Moo = 1*07) flow. The global convergence 


histories for these two cases are also shown in fig- 
ures 4(a) and (b). For illustrative purposes, fig- 
ure 4(a) has two curves; one has all implicit bound- 
ary conditions, and the other has an explicit bound- 
ary treatment on the nozzle surface. The difference 
in the convergence is dramatic and clearly indicates 
the advantage of fully implicit boundary condition 
treatment. 

The transonic case was solved using three differ- 
ent approaches: (1) global iteration only, (2) global 
iteration for the first 160 iterations followed by one 
pass with the hybrid-iteration strategy, and (3) one 
pass through the mesh assuming a fully supersonic 
flow followed by two passes with the hybrid-iteration 
strategy. (The first two passes had a loose (10“ 3 ) 
tolerance, and the final pass had a tight (10~ 12 ) 
tolerance. The results of these three approaches, 
along with the supersonic results, are summarized 
in table 2. 

Table 2. Comparison of CPU Time 
for Dual-Throat Problem 1 

[Tolerance — 10~ 12 ; L 2 ~norm of residual] 


Iteration strategy 

CPU time, sec, for — 

Supersonic 

Transonic 

Global 

68.9 

489 

Global/hybrid 160:1 


381 

Global/hybrid 1:2 


151 

Local 

7.7 



Control Data CYBER 203 time with FORTRAN 
2.1.5 cycle 607 compilation. 


The results again emphasize the fact that local 
iteration is significantly more efficient than global 
iteration for supersonic flows. It is also apparent that 
major gains for transonic flow can be made using the 
hybrid-iteration strategy. 

Another simple problem that indicates the effi- 
ciency of the block-iteration strategy is the Mach 
reflection problem in an inlet-diffuser. The nozzle 
geometry and computational details are described in 
references 14 and 15. A 69 x 31 grid with slight clus- 
tering near the wall and centerline was used. With an 
inflow Mach number of 1.95, a Mach reflection occurs 
(see fig. 5). The flow is almost entirely supersonic, 
with only a small subsonic region existing behind 
the Mach stem. The global convergence of the first- 
order scheme is shown in figure 6. The CPU time 
to obtain a 10“ 12 residual was 48.9 sec, and, with 
the three-pass block-iteration strategy previously de- 
scribed, the CPU time was reduced to 17.8 sec. This 
result, however, was obtained using the solution from 


7 




the previous line as the initial condition, so a further 
improvement could likely be realized. 

A more severe test of the solution strategy is 
described in this case, where the solution of the Euler 
equations in a highly contoured dual-throat inlet is 
sought. The nozzle contour and 69 x 31 grid are 
shown in figure 7. The wall is made up of five 
polynomial arcs. Wall slope is continuous at the arc 
interfaces, and curvature is discontinuous at two of 
the interfaces. 

The flow entering the nozzle is subsonic and orig- 
inates from an infinite reservoir, which is at rest. 
At the inflow boundary, two conditions are obtained 
from evaluation of locally one-dimensional Riemann 
invariants (see ref. 10), which are 

R ± =Vn± — — 7 (23) 

The velocity component V n is the component of ve- 
locity normal to the boundary and is defined to be 
positive in the direction pointing outward from and 
normal to the boundary. For subsonic flow, R~ can 
be prescribed on the boundary and R+ can be eval- 
uated from the interior. The remaining two condi- 
tions that were prescribed at the inflow boundary 
were the entropy (p/p^), which was held constant 
at stagnation conditions, and the y velocity compo- 
nent v , which was set to zero based on a parallel-flow 
assumption at the entrance of the nozzle. The condi- 
tions for the wall and centerline were flow tangency 
and extrapolation of density, temperature, and total 
enthalpy. The flow at the exit station is fully super- 
sonic; thus, all quantities were extrapolated there. 

The initial conditions used to start the calculation 
were relatively simple. The Mach number/area rela- 
tion (assuming quasi-one-dimensional isentropic flow 
and sonic conditions at the first throat) was used to 
obtain the initial entrance conditions. Similarly, ini- 
tial conditions at the exit station were determined 
based on the exit-to-first-throat area ratio flow. The 
Mach number, density, and static pressure were then 
linearly interpolated into the ^-direction. The indi- 
vidual velocity components were chosen such that the 
initial flow direction was parallel to the grid lines in 
the streamwise direction. 

Mach number contours are shown in figure 8. The 
Mach number at the entrance to the nozzle converged 
to M — 0.09 (the Mach number area relations used 
M ^ 0.1). The flow accelerates to sonic conditions at 
the first throat, continues to accelerate downstream 
of the first throat, and a Mach reflection forms in the 
converging region ahead of the second throat. Just 
ahead of the second throat and immediately behind 
the Mach stem, a small subsonic region exists. The 


flow in this region accelerates to sonic conditions at 
the second throat, beyond which the flow becomes 
entirely supersonic. At the exit station, the Mach 
number varied from 1.96 on the centerline to 4.39 on 
the surface. Quasi-one-dimensional theory initialized 
the exit Mach number to 3.22. 

The convergence histories of the first-order and 
fully upwind second-order schemes using global iter- 
ation are shown in figure 9. The sawtooth pattern 
is typical for the higher order upwind methods with 
alternate line Gauss-Seidel sweeping. In the calcu- 
lation, the maximum time step was limited and re- 
sulted in limiting the maximum Courant number to 
405. The results obtained with three hybrid-iteration 
sweeps are presented in table 3. 


Table 3. Comparison of CPU Time for 
Global and Hybrid Iterations 


Difference method 
parameters 

CPU time, sec, for — 

4 > 

% 

Global 

iteration 

Hybrid 

iteration 

0 

N/A 

193 

Not calculated 

1 

-1 

257 

168 

1 

1/3 

259 

171 


The dual-throat problem contains two separate 
subsonic regions which together span 55 percent of 
the x = Constant grid lines. An approximation for 
the optimal performance improvement that could be 
obtained from the hybrid-iteration strategy would 
therefore be 45 percent, since this could be achieved 
only if the supersonic-subsonic interfaces were known 
a priori and if the CPU time associated with the su- 
personic regions were neglected. The results in ta- 
ble 3 show that a 35-percent reduction in CPU time 
was actually obtained. Thus, for this problem, the 
three-pass hybrid-iteration strategy attains approx- 
imately 78 percent of the optimal performance im- 
provement that one could expect with this technique. 

Concluding Remarks 

The family of algorithms presented in this study 
results in an efficient iterative technique for inviscid 
flow problems, particularly in the supersonic regime. 
The new iteration approach for problems contain- 
ing separate regions of supersonic and subsonic flow 
also represents a significant savings in computational 
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effort. The optimal approach to implementing this 
strategy is still unknown. 


NASA Langley Research Center 
Hampton, VA 23665-5225 
September 19, 1986 
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Iteration 

Figure 2. Convergence histories for shock reflection problem. 


Moo =1.10 

NU = 1.07 



Figure 3. Centerline Mach number distribution and wall contour of mildly contoured dual-throat inlet. 




(a) Supersonic; M = 1.10; 69 x 31 mesh. 



(b) Transonic; M = 1.07; 69 x 31 mesh. 

Figure 4. Convergence histories for dual-throat problem. 
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Figure 9. First- and second-order convergence histories for highly contoured inlet problem. 
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